Constraint-preserving Sommerfeld conditions for the harmonic Einstein equations 
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The principle part of Einstein equations in the harmonic gauge consists of a constrained system 
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their accuracy. 
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I. INTRODUCTION 



We present here the implementation and test results of a new formulation of constraint-preserving Sommerfeld 
type boundary conditions for Einstein's equations. The well-posedness of the initial-boundary value problem (IBVP) 
for symmetric hyperbolic systems with maximally dissipative boundary conditions can be established by the energy 
method [l[. An alternative technique, based on the principle of frozen coefficients, Laplace- Fourier decomposition 
and the theory of pseudo-differential operators, can be used to establish well-posedness in a generalized sense even if 
the boundary conditions are not maximally dissipative or the system is not symmetric hyperbolic [2, [3( • This theory 
■ has recently been applied to formulate a well-posed, constraint-preserving IBVP for the harmonic Einstein equations 
with boundary conditions of the Sommerfeld type Q. In this paper, we show how these new boundary conditions 
can be implemented in a finite-difference harmonic code in which the Einstein equations are reduced to second order 
wave equations. The test results presented here show that this new approach has potential value for the computation 
of gravitational waves in a highly dynamical and nonlinear regime. Although our application here is limited to test 
^ | problems, we expect these techniques to further the recent progress in the simulation of black holes by harmonic 
y; evolution 0,Q . 

£h ■ The first well-posed formulation of the IBVP for Einstein's equations was presented in the pioneering work of 
Friedrich and Nagy @, based upon a quite different formulation of the Einstein equations. The underlying pseudo- 
differential theory and how it leads to a constraint-preserving IBVP for the harmonic Einstein equations which is 
well-posed in a generalized sense is described in In that work, the details were presented for the linearized 
Einstein equations but it was explained how the general pseudo-differential theory extends well-posedness to the 
full nonlinear case. Subject to a certain inequality which is necessary to establish the required estimates, there is 
considerable freedom in the detailed form of the Sommerfeld- type boundary conditions. 

In Sec. we describe how this new boundary treatment can be implemented in a fully nonlinear code for the 
simplest choice of the boundary conditions considered in [J]. The resulting scheme is attractive for numerical use. 
In a previous study [l(J, constraint-preserving boundary conditions for the harmonic Einstein equations were based 
upon a combination of Dirichlet and Neumann conditions, which are only marginally dissipative. The description 
of a traveling wave requires the proper inhomogeneous Dirichlet or Neumann boundary data for the wave to pass 
through the boundary. However, in numerical simulations, such inhomogeneous Dirichlet or Neumann data can only 
be prescribed for the signal and the numerical error is reflected by the boundary and accumulates in the grid. Test 
results show that this can lead to poor performance in the simulation of highly dynamical and nonlinear solutions of 
Einstein's equations [lll.ll2j|. For such computational purposes, it is more advantageous to use a Sommerfeld condition 
which is strictly dissipative and allows numerical error to leave the grid [l3| . 

The numerical implementation of the boundary conditions, described in Sec. IIII1 is carried out using two distinct 
approaches. One approach is based upon summation by parts (SBP), which incorporates semi-discrete versions of the 
conservation laws obeyed by the principle part of the system. The stability of the finite difference scheme then follows 
from a discretized energy argument. The second is based upon the embedded boundary method, which is more easily 
applied to the case of a curved boundary. Theorems regarding the stability of the embedded boundary method for 
the second order wave equation have been given for the case of Dirichlet and Neumann boundary conditions fLU [l5| . 
We expect that these theorems can be extended to the strictly dissipative Sommerfeld case. 
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We compare these two numerical implementations using the standardized Apples with Apples linearized wave, gauge 
wave [l6[ and shifted gauge wave [llj tests, modified to include a boundary as documented in [17l| . These tests allow 
a direct comparison of the SBP and embedded boundary methods in a context where the boundary is aligned with 
the grid. The test results are given in Sec. HVl 



II. CONSTRAINT-PRESERVING SOMMERFELD BOUNDARY CONDITIONS 

Our results apply to the generalized harmonic formalism including harmonic gauge source terms [18[ and constraint 
adjustments, as described in [ill. Il9j|. Generalized harmonic coordinates x a = (t,x l ) = (t,x,y,z) are independent 
solutions of the curved space scalar wave equation, 

□z" = -^d a (V=gg a0 d^) = -f", (2.1) 

with gauge source terms T^(x a , g a p) depending on the coordinates and the metric. In terms of the connection T^, 
these harmonic conditions take the form 



where 



C := - r" = 0. (2.2) 



r " = 9 aP K B = ~^d al ^ (2.3) 



and 7 MI/ = yj—gg^"- These C M are the constraints of the harmonic formulation. 
Constraint adjustments of the form 

A itu _ CPA% u (x a ,g af3 , d 7 g af3 ). (2.4) 
can be introduced to modify the reduced system of harmonic equations. These reduced equations then take the form 

:= - VW^ + ^g^X7 a C a + , (2.5) 

which are equivalent to Einstein equations = when the constraints C p are satisfied. When the constraint 
adjustments and gauge source terms vanish, they reduce to the standard harmonic reduction of the Einstein tensor 
(see e.g. 0, H), 

W v ._ GP ,u _ v (^) + l g ^y a T« = 0. (2.6) 

The systems (|2.5I) and (|2.6I) have the same principle part and they both constitute a constrained system of quasilincar 
wave equations with a well-posed Cauchy problem. In terms of the metric, these quasilinear wave equations (12. 5|) can 
be put in the form 



2y^E^ = g alJ d a d 0l ^ + = 0, (2.7) 

where S^ v are terms which do not contribute to the principle part. 

The solutions of the generalized harmonic evolution system (|2.5p are solutions of the Einstein equations provided 
the constraints C M are satisfied. The Bianchi identities, applied to ()2.5() . imply that C M obeys the homogeneous wave 
equation 

V Q V Q C^ + R»C - 2V„(C p An = 0. (2.8) 



The well-posedness of the Cauchy problem for (|2.8|l enforces the unique solution C M = in the domain of dependence 
of the initial Cauchy hypersurface S provided the Cauchy data 7 M "|s and dt^^ v \s satisfy C^\s — dtC^\s — via ()2.5jl . 
It is straightforward to show that these initial conditions are satisfied if the data on S satisfy the Hamiltonian and 
momentum constraints G* = and the initial condition C 1 = 0. 

In order to extend constraint preservation to the IBVP with boundary B it is sufficient to prescribe boundary 
conditions for 7'*" which imply a maximally dissipative homogeneous boundary condition for C M . In [Tol. fl2l|. this was 
achieved by a combination of Dirichlet and Neumann boundary conditions on the various components of 7^" , which 
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then induce a combination of homogeneous Dirichlet and Neumann conditions on the components of C M . However, 
test results showed that Dirichlet and Neumann boundary conditions, which are only borderline dissipative, were 
considerably less accurate than a strictly dissipative Sommerfeld condition [Tl| - These test results were based upon 
exact solutions for which the correct Sommerfeld data were known. Here we consider another approach in which 
constrained Sommerfeld data may be applied consistently in the absence of an exact solution. This Sommerfeld data 
for the components of lead to the homogeneous Dirichlet condition 

C|s = (2.9) 

on the constraints, which is sufficient to guarantee a constraint-preserving well-posed IBVP. It is possible that a 
homogeneous Sommerfeld boundary condition on the constraints would lead to better constraint preservation in 
numerical applications. This requires a second differential order boundary condition on the metric variables for which 
existing theory gives little guidance. See [2l| for a fuller discussion and some promising results. 

Well-posedness depends only on the principle part of the quasilinear system (|2.7p . The pseudo-differential theory 
implies the principle of frozen coefficients by which well-posedness for the nonlinear problem can be established by 
treating the metric g a ^ governing the wave operator in (|2.7| as a constant. Thus the metric can be transformed by 
a linear transformation into Minkowski form and without loss of generality it suffices to establish well-posedness for 
the system of flat space wave equations 

r]^d a d^ v = (2.10) 

subject to the constraints (|2.2[) . In this local Minkowski frame, where \J—g = 1, the constraints reduce to 

:= -d v ^ v - D" = 0. (2.11) 

We choose our local frame so that the +a;-direction is the outward normal to B and the x a = (t, y, z)-directions are 
tangent to B. We write x A = (y, z). The Sommerfeld boundary condition on a scalar field <E> then takes the form 

{d t + d x )<S>\ B = q(x a ), (2.12) 

where q(x a ) represents the prescribed Sommerfeld data. Because the Sommerfeld condition is strictly dissipative it 
leads to a well-posed IBVP in the case of the scalar wave equation. 

Constraint-preserving boundary conditions for the system of linearized Einstein equations (|2.10|) - (|2.1ip can be 
expressed as a hierarchy of Sommerfeld boundary conditions. There are numerous options in this approach Q and 
here we consider the mathematically simplest scheme. 

First we require the 6 Sommerfeld boundary conditions 

(d t +8 x ) 7 ' 4s = ^ i V) (2.13) 
(d t +8 x ) ( 7 tA - 7* A ) = ? ty V) ~ <TV) (2.14) 
(d t + d x ) (j tt -2j tx +j xx ) = q u {x a )-2q tx {x a ) + q xx (x a ), (2.15) 

where the g's are freely prescribed Sommerfeld data. Next, the constraints are used to supply 4 additional boundary 
conditions in the hierarchical order 

C A \ B = -d a At - d xl Ax - d Bl AB - f A (x a ) = (2.16) 
€% -C X \ B = -d t {l u - l xt ) - d x (j xt - 7 XX ) - d B ( 7 tB - T * B ) - f l {x a ) + f x (x a ) = (2.17) 
C% =-da u -d xl tx -d B l tB -f t (x a )=0. (2.18) 

By using (|2.13[) - (|2.15[) . the boundary conditions (|2.16|) - (|2.18[) can be re-expressed in the Sommerfeld form 

C A \ B = -\{d t + d x ){ 1 At + 1 Ax ) - d t (j At - J Ax ) - d Bl AB + l -{q tA - q xA ) -f A = Q (2.19) 

C%-C X \ B = -±(d t + d x )( 1 tt - 1 xx )-d t ( 1 tt -2 1 xt + 1 xx ) 

~ d B (~f tB - 7 xS ) + i((? t4 - 2q tx + q xx ) - f 4 + t x = (2.20) 
C% - -~(d t + 9 x )( 7 t4 + l xx ) - a t ( 7 44 - l tx ) + \{q tt - 2q tx + q xx ) - d Bl tB - f 4 = 0. (2.21) 
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Equations (|2.13[1 - (|2.15|) and (|2.19|) - (|2.21[) form a hierarchical sequence of inhomogeneous Sommerfeld boundary condi- 
tions in which the source terms for (|2.19[) - (|2.2ip are provided by previous members in the hierarchy. In the linearized 
problem, the boundary conditions (|2.13p - (|2.15|) . along with compatible initial data that satisfy the constraints, deter- 
mine unique solutions of the wave equations for j AB , j At — j Ax and 7" — 2-f tx + j xx . These then provide the source 
terms in (|2 . 1 9[) and (|2.20[) , which determine unique solutions for r y At + ^ Ax and 7" — "f xx . Finally, the source terms 
in (I2.20|) can be determined to provide a unique solution for 7" + j xx . 

The ten linearly independent Sommerfeld-type boundary conditions (|2.13|) - (|2.15p and (|2.19p ~ (|2.2ip for the compo- 
nents of 7^ give rise to a unique solution of the constrained linearized problem. However, it is important to emphasize 
the following points: 

1. Proof of the well-posedness of the IB VP requires estimates on the derivatives of the solution at the boundary. For 
that purpose, it is required to apply the pseudo-differential theory to construct a symmetrizer in Fourier-Laplace 
space, as described in Q. 

2. The pseudo-differential theory allows well-posedness to be extended locally in time to the nonlinear IBVP for 
the harmonic Einstein's equations, where the relation 7^ = ^/—gg tJ/l> converts (|2.7p into a quasi-linear equation. 

3. Sommerfeld boundary conditions are strictly dissipative but they completely remove all reflections only in highly 
idealized cases. 

4. The hierarchical structure of the boundary conditions, together with their dissipative property, is a very good 
heuristic procedure to formulate a stable finite difference approximation but by itself does not provide a proof 
of stability. 

5. In the electromagnetic case, our approach gives rise to constraint-preserving Sommerfeld-type boundary condi- 
tions on the vector potential which can be given a physical interpretation in terms of the Poynting vector Q • An 
analogous interpretation does not exist in the gravitational case. When q AB has vanishing trace-free part, the 
boundary condition (|2.13p implies that the outgoing null hypersurfaces emanating from the boundary have van- 
ishing shear. This is related to the vanishing of incoming radiation but in a very gauge dependent way. Friedrich 
and Nagy |9( stress this same caveat concerning their boundary conditions, which specify the Newman-Penrose 
Weyl component Wo associated with the outgoing null hypersurfaces. 

6. In 0], well-posedness of the IBVP was established for the special case where the boundary x — contains the 
timelike normal to the initial Cauchy hypersurface t — 0. More generally, the initial Cauchy hypersurface is 
given by T = t — kx so that the initial data is different from the initial data at t = 0. Well-posedness of this 
more general IBVP can be reduced to well-posedness of the special case by first considering a pure Cauchy 
evolution from T = to t = 0. The general case corresponds to a boundary that is a "moving" with respect to 
the Cauchy hypersurface. This introduces a shift term in the wave operator. In the next section, we give the 
details of how to translate the conditions (|2.13p - (|2.15p and (|2.19p - (|2.2ip into a nonlinear computational scheme 
for a general metric which is not in the shift-free Minkowski form. 



III. NUMERICAL IMPLEMENTATION 

Our implementation is based upon the Abigel code [l(| E3] , m which the harmonic system (|2.7I ) is integrated in 
the first order in time form 

dtT ^ = i^( 7)9 . 7)Tifil . T ). (3.2) 

Here x l = (x, y, z) and <9; denotes spatial derivatives. The code is based upon an explicit, second order accurate 
finite-difference scheme. Introduction of a spatial grid and finite difference approximations for the spatial derivatives 
reduces (I3.ip ~ (|3.2p to a large set of ordinary differential equations (method of lines), which are evolved with a fourth 
order Runge-Kutta integration. Details of the finite difference approximations are given in where the system 

(I3.ip - (|3.2p was expressed in flux conservative form and summation by parts (SBP) was used to apply the boundary 
conditions, which enforced the semi-discrete version of the conservation laws obeyed by the principle part The tests 
in this paper have been carried out with the W form of the algorithm described in [12J, in which certain nonlinear 
coefficients are approximated by their averages between grid points. We also add artificial dissipation to (|3.I p and 
p.2p by the modifications 

W dff + (3-3) 
d t T» v -> d t T^ + e T V A T^ v , (3.4) 
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where T> 2 is the SBP approximation for the Laplace operator. 

The SBP approach requires that the boundary be aligned with the numerical grid. Here we also consider imple- 
mentation of the boundary conditions by means of the embedded boundary method [3, El , which is applicable even 
when the boundary is not aligned with the grid. 



A. The Sommerfeld boundary data 



We describe the implementation of constraint-preserving Sommerfeld boundary conditions for the fully nonlinear 
harmonic system (|2.5j) . For purpose of discussion, we locate the boundary B at x = const with the outer normal in 
the +ir-direction. We denote the coordinates intrinsic to the boundary by x a = (t,y,z) = (t,x A ). We introduce an 
orthonormal tetrad (T M , X^ 1 , Y^, Z<- L ), oriented at B so that X^ = (0, l/y/g xx ,0, 0) is the unit outward normal and 
= (1/^—gtt, 0, 0, 0) is a unit vector in the evolution direction. In addition to the standard requirements of the 
IBVP that the Cauchy hypersurfaces be spacelike and that the boundary be timelike, we require that T M be timelike, 
i.e. that the evolution be subluminal. 

The Sommerfeld data for 7^ consist of 

K a d a ^ v = q» v , (3.5) 

where = T p + X^ and K^K p — 0. The description of the free and constrained components of is algebraically 
more complicated than for the Minkowski case in Sec. [Til For this reason, it is useful to introduce the projection 
operators 

M; =51 + T p T v - X p X v = Y p Y v + Z p Z v , (3.6) 
which projects vectors into the (Y~ M , Z^) spatial-plane tangent to the boundary, and 



K=K + n^K", (3.7) 



where = — X^ with K^L p = —2, which projects vectors into the null 3-space spanned by (L M , F M , Z^). We 
have P p K^ = and P p L v = 0, i.e. P p projects forms W p into the subspace orthogonal to K p . We have 

P; = Ml - \k^L v . (3.8) 

The freely prescribed Sommerfeld data corresponding to (|2.13l) - (|2.15[> consist of the projection 

QV, = p^pVcfP. (3.9) 

The constrained data consist of = P£Lpq al3 , analogous to (|2~T9| - ([2~20| : and Q = L Q i /3 q Q/3 , analogous to (j2~2Tj) . 
In order to determine the constrained data, we use the identity 

d a T u = -\L a q a » + P p a d pl av - (3.10) 



The constraints then imply 



so that 
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l^f" = -\L a iT + P p ^ P l av (3.11) 



= 2P£pPd pl afi + 2^P^T a (3.12) 



and 



Q = 2L a P p p d pl af} + 2^L a t a . (3.13) 
The full set of Sommerfeld data consist of 



5' 



1 



qh» _ q(p-k u) + -QK^K". (3.14) 
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As in (|2.19[) - (|2.21[) , the derivatives normal to the boundary occurring in (|3.12[) and (|3 . 1 3|) can be eliminated via the 
identities 

PZPpdpl^ = X a Q» a - P^K p T p d pl a ^ + P p M p d pl a0 (3.15) 

and 

L a P p d pl ap = \K a Q a - L a K p T p d pl ap + L a M p d pl a0 . (3.16) 

Combined with (|3.12[) and (|3 . 1 3|) . these equations determine the constrained boundary data Q p and Q in terms of 
quantities intrinsic to the boundary, 

= 2^K a Q pa - P£K p T p d pl af3 + P£M p d pl afi + V^gPgfA (3.17) 

Q = 2 {^K a Q a - L a K p T p d pl a(} + L a M p d pl afJ + ^L a f a ^j . (3.18) 

B. The boundary update algorithm 

In Sec. MI Al we have described how the boundary values of 7 A " / , T^ v — dt^ v and the free Sommerfeld data Q^ v 
at time t determine the full Sommerfeld data at time t. We now consider how these quantities are updated. The 
update of the boundary values of ^ v to the next time step t + At is straightforward via ()3.ip ). The update of is 
more complicated. Here we use two distinct algorithms, the SBP algorithm and the embedded boundary algorithm, 
for updating the boundary values of T Miy . This provides two competitive update algorithms, which are compared in 
test problems in Sec. IIVI 

SBP algorithm. This is the algorithm described in [12]. The evolution equation (|3.2p is applied at the boundary 
points to update T M ", with field values at the resulting ghost points eliminated by the boundary condition in a manner 
that enforces discrete conservation laws. Full details for the cases of Dirichlet, Neumann and Sommerfeld boundary 



conditions are given in 12j. The only new ingredient here is the use of constraint-preserving Sommerfeld conditions. 
However, because the constrained Sommerfeld data and Q can be updated after the update of r y p " and ', there 
is no essential change in the numerical algorithm. 
Embedded boundary algorithm. 

The SBP algorithm uses the fact that the boundary is aligned with the grid, so that boundary points are grid 
points. In the case of a spherical boundary and a Cartesian grid, this is generally not the case unless multi-block 
techniques are used. An alternative approach for treating a curved boundary with a Cartesian grid is the embedded 
boundary method, in which field values at ghost points are updated from interpolations using the boundary data, as 
opposed to applying the evolution equation. 

Assume that we are given the values of and T^ v at time t and the free data Q^ u at all times. As above, this 
determines the full constrained Sommerfeld data q pu at time t. After updating the boundary values of 7^ to time 
t + At, we update the boundary values T^ v using 

= K a d a ^ v = K l T^ + K^n^, (3.19) 

where the spatial derivatives K t dij )J " 1 ' at the boundary are determined to by an interpolation scheme (see below). 
When all components of the Sommerfeld data q pu are supplied analytically, this provides updated boundary values 
for , which we refer to as the analytic version of the embedded boundary algorithm. 

In the absence of an analytic solution, only the unconstrained components of the Sommerfeld data can be 
prescribed freely. The update of is then first carried out for those components determined by the free data, i.e. 

r pv = P^P v p T a0 7 (3.20) 

for which (|3 . 1 9[) gives 



Q pv = K 1 ^ + P^PlK i d ll afi . (3.21) 
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In order to update the boundary values of the remaining components of T^", we must first update the constrained 
Sommerfeld data and Q by expressing all time derivatives on the right hand sides of (|3.17[) and (|3.18p in terms of 
previously updated quantities. From (|3.17p . we obtain 

Q» = 2 (\K a Q^ a - —r== K a T^ a + M l a T^ + P£Mld tl aP + y/^P^fA . (3.22) 

This determines 



t» = P£L fj T ai:i (3.23) 

through 

_ K t T n + p^ L/3 K t d l -/ a(3 . (3.24) 



Similarly from (|3.18p . we next obtain 



9tt 



K a r a + M l a T a + L a Mldii 



a/3 



-gL a f a 



(3.25) 



This determines 



t = L a L p T al3 (3.26) 

through 

Q = K 1 t + L a L p K l da a(3 . (3.27) 



These pieces allow us to construct 

T ^ = T ^ _ T {n K v) + \ T K»K V . (3.28) 

In carrying out this update of T^ v it is necessary to compute the boundary values of the spatial derivatives K l dij ,J ' 1 ' 
which appear in (|3.19|) . This is accomplished by an interpolation scheme patterned after the embedded boundary 
treatment for a Neumann condition given in [l5|| . In this scheme, the value of K l di^ fil ' at a boundary point B, with 
grid values x l B = (Ib, Jb, Ks)h, is obtained from the 1-dimensional Lagrange polynomial for determined by B 
and M interior points, C m (1 < m < M), lying along the curve 

x i =x i B + —(x-x B ) (3.29) 

normal to the boundary. We choose (x B — xc m ) = mh so that the points C m lie in (y, z) planes through the a;-grid. The 
values of at the points C m are obtained by 2D interpolations in those planes. We use a 2D Lagrange interpolant 
based upon a stencil of grid points which is symmetrical about the point (i/b,zb)- In order to avoid extrapolation, 
the size of the stencil has to be adjusted to the size of K l /K x . For all 2D test cases considered in Sec lIVl it suffices 
to use a 3 x 3 (or larger) stencil in the planes determined by C m . 

Note that this scheme differs from the straightforward computation of K l di^^ v by centered differencing in the y 
and z directions and one-sided differencing in the x direction, which would give a qualitatively different and less 
accurate approximation. However, in the simple case where K y — K z — 0, this interpolation scheme does reduce to 
approximating K l di r ~i p,v — K x d x ^ v by a one-sided finite difference at B. For a second order accurate scheme, only 
3 points C m (M — 3) are necessary to construct the Lagrange polynomial approximating (|3.29[) . However, with a 
second order accurate interior evolution algorithm, the error would then be largest at the boundary because of the 
use of a 1-sided derivative. On the other hand, it is simple to use more points if higher accuracy at the boundary 
is desired. In this respect, the embedded boundary algorithm is more flexible than the SBP algorithm, for which 
the order of accuracy of the boundary algorithm is coupled with the order of accuracy of the interior algorithm. In 
the embedded case, for any internal accuracy, the order of accuracy of the boundary condition can be made high as 
desired at little additional computational expense. In the tests results shown in the next section, we use M = 5. We 
also add a corrector step to the boundary update algorithm to improve accuracy. 
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IV. TESTS 



We conduct three boundary tests which have been proposed for the AppleswithApples ( AwA) test suite [Tj] . These 
tests extend the original linearized wave, gauge wave and shifted gauge wave tests with periodic boundaries, i.e. a 
3-torus T 3 without boundary, to include non-trivial boundaries by opening up the :r-axis of the 3-torus to form a 
manifold with smooth T 2 boundaries at x = ±.5. In this way, the tests avoid the complication of sharp boundary 
points. The corresponding metrics are 



Linearized wave: 



Gauge wave: 



• Shifted gauge wave: 



where in all cases 



ds 2 = -dt 2 + dx 2 + (1 + H)dy 2 + (1 - H)dz 2 , (4.1) 
ds 2 = (1 - H)(-dt 2 + dx 2 ) + dy 2 + dz 2 , (4.2) 
ds 2 = -dt 2 + dx 2 + dy 2 + dz 2 + Hk a k fs dx a dxP, (4.3) 



H = H(x-t) = Asm( 27r( ^ t} ) , (4.4) 



and 



k a =d a (x-t) = (-1,1,0,0). (4.5) 



These metrics describes sinusoidal traveling waves of amplitude A propagating along the x-axis. Two dimensional 
features are tested by rotating the coordinates according to 

x = -^=(x -y), y = -j={x' + y'). (4.6) 

which produces a wave propagating along the diagonal. 

The linearized wave test is run with an amplitude A = 10~ 8 . It is most efficient for revealing problems arising from 
nonlinearity to run the gauge wave and shifted gauge wave tests with amplitude A — 0.5. In some cases we also run 
with smaller amplitudes in the range A = .01 to A = .1 (the original AwA specifications) to reveal the emergence of 
nonlinear features. In all other respects, we retain the original AwA specifications: 

• Wavelength: d = 1 in the ID simulation and d! = 1/V2 in the 2D simulation. 

• Simulation domain: 

ID: x E [-0.5, +0.5], y = 0, z = 0, d=l 

diagonal: x <E [-0.5, +0.5], y G [-0.5, +0.5], z = 0, d' = y/2 

• Grid: x n — —0.5 + ndx, n = 0, 1 . . . 50,o, dx = dy = dz = l/(50p), p = 1, 2,4 

• Time step: dt = dx/4 = 0.005/p . 

The grids have N = 50p = (50, 100, 200) zones. (At least 50 zones are required to lead to reasonable simulations for 
more than 10 crossing times.) The ID tests are carried out for t — 1000 crossing times, i.e. 2 x 10 5 p time steps, and 
the 2D tests for 100 crossing times. 

As an example of how the Sommcrfcld boundary data is prescribed, consider the ID shifted gauge wave. The 
Sommerfeld operator at the right boundary x = +.5 is given by K a d a = vl — H (dt + d x ) and the corresponding 
Sommerfeld data vanishes (q^ = 0). At the left boundary x = — .5, where the outward normal is in the (— a;)-direction, 
the Sommcrfcld operator is 

K a d a = (T a - X a )8 a = \ +H d t - VT=Hd x (4.7) 

V 1 ~ ti 
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and the resulting non-vanishing components of Sommerfeld data are 

q tt = q t x=qXX = __^_ dtH{x _ t) _ (4g) 

With this inhomogeneous Sommerfeld data, the wave enters through the boundary at x = — .5, propagates across the 
grid and exits through the boundary at x — +.5. 

At the left boundary x — — .5, the inward null direction — + X^, which enters the projection operator (|3.7p . 
satisfies L^q^" = 0. Consequently, the non- vanishing components of the free data Q^ v are also given by (|4. 8|) and 
the constrained components vanish, i.e. = Q = 0. Of course, in a simulation using the constraint preserving 
Sommerfeld algorithm, numerical error gives rise to non- vanishing values for and Q. Thus it makes a difference 
at the numerical level whether Q M and Q are given their analytic values (zero) or their values produced by the 
constraint-preserving algorithm. These alternatives will be compared in carrying out the tests to provide a measure 
of the efficacy of the constraint-preserving algorithm. 

Thus for each of the two boundary algorithms (SBP and embedded boundary) we run the tests (i) with all f com- 
ponents of Sommerfeld data q^ v provided by the analytic solution and (ii) with the 6 free components of Sommerfeld 
data Q^ v provided by the analytic solution and the remaining components and Q provided by the constraint- 
preserving algorithm. Wc distinguish the corresponding tests results by the labels ASBP and AEMB, respectively, for 
the SBP and embedded algorithms with fully analytic data; and CSBP and CEMB for the corresponding tests with 
constraint-preserving data. 

In the ID gauge wave tests, the curves (|3.29l) normal to the boundary pass through the grid points on the x-axis, 
so that the interpolations required for the embedded boundary algorithm are trivial. The 2D gauge wave and shifted 
gauge wave both provide a non-trivial test of the interpolation scheme. For the case of the shifted gauge wave, the 
curves (|3.29[) at the boundaries xb = ±-5 are given by 

K y 

V = Vb + jt^(x - x B )- (4.9) 

Harmonic gauge forcing terms were not found to be effective in the boundary-free gauge wave tests and we have not 
included them in the present tests. (Gauge forcing is important in spacetimes where harmonic coordinates become 
pathological, e.g. the standard i-coordinate in Schwarzschild spacetime is harmonic but singular at the horizon.) 
Constraint adjustments were not necessary except in trying to stabilize the constrained SBP shifted gauge wave runs. 

We use the £oo norm to measure the error 

£(4) = ||$p-$ana||ao (4-10) 

in a grid function <& p with known analytic value $ Q „ a - We measure the convergence rate at time t 

using the p = 2 and p — 4 grids (N = 100 and N = 200). (The p = 1 grid is not very accurate and is mainly used for 
debugging.) For convergence studies it is also useful to graph the rescaled error 

£p = JgP p -$a«a||oc, (4.12) 

which is normalized to the p — 4 grid. 



A. Linearized wave tests 



Table U shows the convergence rates of the error in g yy for the ID linearized wave, measured at t = 10 and t = 50 
crossing times. Clean second order convergence is maintained for all four algorithms, irrespective of whether the 
complete Sommerfeld data is supplied from the analytic solution (ASBP and AEMB), or whether it is constrained 
(CSBP and CEMB). At 1000 crossing times, the four algorithms continue to give excellent agreement with the analytic 
solution. The graphs in Fig. [1] show excellent phase agreement and a small difference in amplitude at t = 1000 in the 
comparison between the analytic solution for g yy (x) and the results the CEMB and CSBP results. On the scale of 
Fig. [TJ the graphs for the AEMB and ASBP algorithms are indistinguishable from the analytic solution. 

Table |TT] shows the convergence rates of the error in g yy for the 2D linearized wave test, measured at t = 1, t = 5 
and t = 10. Second order convergence is cleanly maintained for the AEMB and ASBP algorithms. The convergence 
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ALGORITHM 


t=10 


t=50 


AEMB 


2.04 


2.05 


ASBP 


1.99 


1.93 


CEMB 


2.04 


1.97 


CSBP 


1.99 


2.03 



TABLE I: Convergence rates of £{g yy ) for the ID linearized wave test, amplitude A = 10 




-0.4 -0.2 0.2 0.4 0.6 

X 



FIG. 1: The graphs compare the performance of the CEMB and CSBP algorithms for the ID linearized wave test on the p = 4 
grid with the analytic solution. The CEMB and CSBP snapshots of g yy {x) — 1, shown at t — 1000, are indistinguishable and 
show only a small amplitude discrepancy with the analytic solution. 



rates of the CEMB and CSBP show some deterioration from second order at t = 10 as the truncation error from the 
boundary algorithms accumulates. Nevertheless, at t = 100, the 2D results for all the algorithms remain in excellent 
agreement with the analytic solution. The projection operators used in the constrained CSBP and CEMB algorithms 
introduce the small errors shown in Fig. [2j where snapshots of g yy — 1 are graphed for y = and t = 100, using the 
p — 4 grid. Overall, the linearized wave tests show that both the SBP and embedded algorithms give excellent results 
for either fully analytic or constrained Sommerfeld boundary conditions. 



ALGORITHM 


t=l 


t=5 


t=10 


AEMB 


1.98 


2.00 


1.97 


ASBP 


1.96 


1.93 


1.92 


CEMB 


1.98 


1.94 


1.62 


CSBP 


1.94 


1.89 


1.77 



TABLE II: Convergence rates of £(g yy ) for the 2D linearized wave test. 



B. Gauge wave tests 

For the ID gauge wave tests, Table Hill shows the convergence rates of the SBP and embedded algorithms, measured 
at t = 10 and t — 50 crossing times. Second order convergence at t = 10 is clean in both cases. At t = 50, the 
convergence rates of the embedded algorithms show slight deterioration, for the reasons explained below. 

All the algorithms remain stable for 1000 crossing times for the ID runs on the p = 4 grid, without use of constraint 
adjustments or artificial dissipation except for a small amount of dissipation, (|3.4[) with = 0.01, for the constrained 
CSBP algorithm. Figure [3] shows the long term performance of the ASBP and AEMB algorithms. Both maintain 
excellent accuracy for 1000 crossing times. For the ASBP algorithm there is negligible error growth. For the embedded 
AEMB algorithm, there is long wavelength error corresponding to the harmonic instability [l3| 



ds\ = e At (l - H)(~dt 2 + dx 2 ) + dy 2 + dz 2 . 



(4.13) 
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FIG. 2: Snapshots of g yy {x) — 1 for the 2D linear wave tests obtained at t = 100, setting y = 0. On the left plot, the snapshots 
of 

9yy ~ 1 f° r the CEMB and CSBP algorithms are compared with the analytic solution. On this scale the errors are almost 
imperceptible and a zoom is given in the the right plot. 

of the gauge wave spacetime (|4.2[) . For any value of A, (|4.13p is a flat metric which obeys the harmonic constraints. 
As depicted in Fig. [3l the resulting profile of the AEMB error contains two peaks. The positive peak dominates in the 
beginning, but there is an overall drop in the waveform, due to the instability, which causes the error to pass through 
a minimum and then grow again as the negative peak dominates. The SBP algorithm is designed to suppress this 
instability by means of discrete conservation laws for the principle part of the evolution equations. The embedded 
algorithm excites the instability although at a fairly innocuous level. 

Similar ID results hold for the constrained algorithms CSBP and CEMB, as shown in Fig.[U The use of a 5 point 
(M = 5) Lagrange polynomial is essential for the good performance of the constrained algorithms. Figure [5] shows 
the rapid error growth which would result from the use of 3 or 4 points, again arising from excitation of the long 
wavelength instability (|4. 13|) . 



ALGORITHM 


t=10 


t=50 


AEMB 


1.97 


1.87 


ASBP 


2.00 


2.00 


CEMB 


1.97 


1.87 


CSBP 


1.98 


1.95 



TABLE III: Convergence rates of £(g xx ) for the ID gauge wave test. 




100 200 300 400 , 500 600 700 800 900 1000 -0.4 -0.2 0.2 0.4 

t X 



FIG. 3: The performance of the analytic SBP and embedded algorithms for the ID gauge wave test. On the left, the error 
norm £{g xx ) is plotted vs t. On the right, snapshots of the error in g xx for the AEMB algorithm are shown at t = 10 and 
t = 1000. This error is long wavelength, consisting of two peaks. The positive peak dominates in the beginning, but an overall 
drop in the profile causes the error to pass through a minimum and then to grow again as the negative peak dominates. This 
behavior is due to the long wavelength instability in the gauge wave spacetime. 
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100 200 300 400 ,500 600 700 800 900 1000 

FIG. 4: Plot of the error norm £{g X x) vs t for the ID gauge wave simulations with the constrained SBP and embedded 
algorithms, obtained with the p = 4 grid. 




200 400 , 60(1 800 1000 -0.4 -0.2 0.2 0.4 
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FIG. 5: On the left, the error norm £(g X x) is plotted vs t for the ID gauge wave test using the CEMB algorithm with Lagrange 
polynomials based upon M — 3, 4 and 5 points. Only for M — 5 is there no long term error growth. On the right, the snapshots 
of the error in g xx at t = 1000 show the long wavelength mode which is excited in the M = 3 and 4 cases. 

The convergence rates for the 2D gauge wave tests shown in table IIVI indicate clean second order convergence 
up to t = 10. The graphs of the error in Fig. [5] show that the analytic SBP and embedded algorithms maintain 
excellent accuracy up to t = 100. However, the constrained algorithms excite the long wavelength instability (|4.13|) 
at t ~ 55 for CSBP and t » 60 for CEMB. Neither numerical dissipation nor constraint adjustment lead to significant 
improvement. Higher order accuracy of the boundary condition also does not seem to help. 



ALGORITHM 


t=l 


t=5 


t=10 


AEMB 


2.02 


2.03 


2.03 


ASBP 


2.02 


2.03 


2.02 


CEMB 


2.02 


2.01 


2.02 


CSBP 


2.02 


2.01 


2.03 



TABLE IV: Convergence rates of £(g xx ) for the 2D gauge wave test. 



C. Shifted gauge wave tests 

For the ID shifted gauge wave tests, Table [V] shows the convergence rates of the SBP and embedded algorithms 
measured at t = 10 and t = 50 crossing times. Second order convergence is fairly clean at t — 10, with some drifting 
of the convergence rates evident at t = 50. 
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FIG. 6: Plots of £(g X x) vs t for the 2D gauge wave tests. The left graphs compare the error for the analytic SBP and embedded 
algorithms. The right graphs, which compare the error in the constrained SBP and embedded algorithms, exhibit the excitation 
of the long wavelength instability. 



ALGORITHM 


t=10 


t=50 


AEMB 


1.97 


1.89 


ASBP 


2.08 


2.34 


CEMB 


1.90 


1.97 


CSBP 


2.08 


2.25 



TABLE V: Convergence rates of £(g X x) for the ID shifted gauge wave test. 



A long wavelength instability also exists in the shifted gauge wave spacetime (|4.3[) [Tlj . 

ds\ = -dt 2 + dx 2 + dy 2 + dz 2 + (h - 1 + e A *^ k a k l3 dx a dx f} , (4.14) 

where 

Ad (2ir(x-t)\ 

Although this metric does not solve Einstein's equations, for any value of A it satisfies the standard harmonic form (|2.6p 
of the reduced Einstein equations, i.e. the equations governing numerical evolution without constraint adjustment. 
This instability is the major source of error. Figure [7] exhibits the long term performance for runs with the analytic 
SBP and embedded algorithms. They maintain excellent accuracy for 1000 crossing times, although the snapshots 
show that the embedded algorithm has produced a low level excitation of the instability. 

The larger boundary errors of the constrained algorithms leads to stronger excitation of the long wavelength 
instability. Without constraint adjustment, the runs on the p — 4 grid crash at about t rs 88 for the CSBP algorithm 
and t w 103 for the CEMB algorithm. Numerical dissipation was necessary in the CEMB run but had insignificant 
effect on controlling the long wavelength instability in either of the constrained algorithms. However, constraint 
adjustment was moderately effective. Of the various adjustments of the form (|2 .4[) which were considered in [ll[, the 
longest constrained runs were obtained with the choice 

bC a V t 

Altv = b > (4 16) 



Here 



e pa =g P a~ ^(V p t)V ff i (4.17) 

is the natural metric of signature (+ + H — h) associated with the Cauchy slicing. This adjustment with 6=1 extended 
the CSBP run to t w 155 and the CEMB run to t w 183 crossing times, as indicated in Figure [HI The response to 
constraint adjustment is evidence of the constraint violating origin of the analytic instability. The instability is also of 



14 




FIG. 7: Performance of the analytic SBP and embedded algorithms for the ID shifted gauge wave. On the left, the error norm 
£(g xx ) is plotted vs t. On the right, snapshots of the error for the AEMB algorithm with the p — 4 grid are shown at t = 10 
and t = 1000. The error is long wavelength, consisting of two peaks. The overall rise in the profile is due to low level excitation 
of the long wavelength instability (I4.14p . 

nonlinear origin, which can be seen from the comparison with the runs of lower amplitude shown in Fig. [5] With lower 
amplitude, the discrete conservation laws of the SBP algorithm begin to control the instability and give performance 
comparable to the CEMB algorithm. As was previously found for shifted gauge wave tests with periodic boundary 
conditions [ll| , constraint damping [ll| introduces oscillations with unacceptably large error and does not appreciably 
suppress the instability. Nevertheless, the performance of the constrained Sommerfeld algorithms is vastly superior 
to the Neumann-Dirichlet constrained algorithm considered in as can be seen in figure [5] 




,1: I r r -r I , I , I , l_ 

20 40 60 80 100 120 140 160 180 



FIG. 8: Plot of £(g X x) vs t for the ID gauge wave with shift simulations for the constrained CEMB and CSBP algorithms, with 
numerical dissipation er = 0.001 and constraint adjustment with 6 = 1, and for the Neumann-Dirichlet constrained preserving 
algorithm. The comparison shows that both Sommerfeld constrained algorithms clearly outperform the Neumann-Dirichlet 
constrained algorithm. 

The convergence rates for the 2D shifted gauge wave tests shown in table lVII indicate clean second order convergence 
up to t = 10. Figure [10] exhibits the long term performance for runs with the analytic and constrained SBP and 
embedded algorithms. For the analytic algorithms (left plot), both the ASBP and the AEMB algorithms maintain 
excellent accuracy for 1000 crossing times. This is one of the few cases where the error in the embedded algorithm is 
larger than the error in the SBP algorithm. This results from the 2D interpolation error introduced by the embedded 
algorithm. For the constrained algorithms (right plot), the long wavelength instability is excited earlier by the CSBP 
algorithm, which crashes at t ss 22, while the CEMB algorithm runs up to t i=a 35. Efforts to control the instability 
by constraint adjustment and numerical dissipation had insignificant effect in prolonging the runs. 
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FIG. 9: Plots of £(g X x) vs t for lower amplitude simulations of the ID gauge wave with shift simulations for the constrained 
versions of the SBP and embedded algorithms. The nonlinear nature of the stability is evident. 



ALGORITHM 


t=l 


t=5 


t=10 


AEMB 


2.00 


2.01 


2.03 


ASBP 


2.00 


2.02 


2.05 


CEMB 


2.00 


2.00 


1.97 


CSBP 


2.00 


2.00 


2.04 



TABLE VI: Convergence rate of £(g X x) for the 2D shifted gauge wave test. 



V. CONCLUSION 



The preceding tests involved linearized wave, gauge wave and shifted gauge wave metrics for which the exact 
solutions provide the correct boundary data. The results provide some definitive conclusions. First, the results show 
that constraint-preserving Sommcrfeld boundary conditions give good long term accuracy, which is vastly superior to 
previous shifted gauge wave test results for constraint-preserving Dirichlet-Neumann conditions. 

Of foremost importance, the analytic versions of both the SBP and embedded algorithms clearly outperformed 
the constrained versions. The knowledge of the exact Sommerfeld boundary data avoids the need for computing 
constraint-preserving boundary conditions. One might have expected numerical noise to generate constraint violating 
error which the analytic algorithms could not properly handle. While this undoubtedly occurs, the analytic algorithms 
nevertheless outperform the constrained algorithms because the prime source of error is due to the long wavelength 
instabilities inherent in the nonlinear test problems. The additional error introduced by computing constrained 
Sommcrfeld data leads to earlier excitation of these long wavelength instabilities. 

Also of importance, our tests results show no clear advantage of either the SBP or embedded treatments of the 
boundary condition. In the more complicated case of a curved boundary, both approaches are being developed by 




FIG. 10: Plots of £(g X x) vs t for the 2D gauge wave with shift. The left graphs compare the analytic SBP and embedded 
algorithms and the right graphs compare the constrained algorithms. 
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their advocates in the computational mathematics community. 

In a realistic problem, such as the binary black hole problem, a global scheme is necessary to provide physically 
correct outer boundary data, either by using hyperboloidal time slices to extend the Cauchy evolution to infinity 
(see [22|, HH for reviews) or by matching to an exterior characteristic or perturbative solution (see [24| for a review). 
Harmonic evolution offers important advantages for Cauchy-characteristic matching (CCM) which have led to suc- 
cessful matching in the linearized regime [Iol | . Primarily, the harmonic constraints can be enforced by propagating the 
Cauchy coordinates on the characteristic grid. This provides the proper Jacobian for injecting the Cauchy boundary 
data from the characteristic solution. The results of this paper should supply helpful guidance for extending the 
harmonic version of CCM to the nonlinear problem. 
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